C
C     adjust the geometry according to the incident wave number
C     and the incident angle
C     if alph, beta, phai and thet are all 0, then it doing nothing
C     
      subroutine geometry_init(nL, r0, k, alph, beta, phai, thet)
      implicit double precision (a-h,o-z)
      include 'gmm01f.par'
      parameter (nmp=np*(np+2),nmp0=(np+1)*(np+4)/2)
C
      double precision k, r0(6,nLp)
      double precision beta, phai, thet

C
C     根据入射角度变换坐标(程序假设平面波x-极化，z-传播!), 重新计算金属球的几何信息
C
      pih = dacos(0.d0)

C     其中alph表示传播方向与+z轴夹角
      alph = alph*pih/90.d0
      ca = dcos(alph)
      sa = dsin(alph)

C     beta表示极化方向与+x轴夹角
      beta = beta*pih/90.d0
      cb = dcos(beta)
      sb = dsin(beta)
      do i = 1,nL
         x0 = r0(1,i)
         y0 = r0(2,i)
         r0(1,i) = cb*x0 - sb*y0
         r0(2,i) = sb*x0 + cb*y0
      enddo

C     convert the angle to 0 ~ pi
      phai = phai*pih/90.d0
      cb = dcos(phai)
      sb = dsin(phai)
      thet = thet*pih/90.d0
      cz  =dcos(thet)
      sz = dsin(thet)
      do i=1,nL
         x0 = r0(1,i)
         y0 = r0(2,i)
         z0 = r0(3,i)
         r0(1,i) = ca*cz*x0 - (ca*sz*sb + sa*cb)*y0
     +        + (ca*sz*cb - sa*sb)*z0
         r0(2,i) = sa*cz*x0 - (sa*sz*sb - ca*cb)*y0
     +        + (sa*sz*cb + ca*sb)*z0
         r0(3,i) = - sz*x0 - cz*sb*y0 + cz*cb*z0
      enddo

      write(6,*) 'Original sperhere poisition:'
      do i = 1,nL
         write(6,'(i5,4f14.5)') i, r0(1,i), r0(2,i), r0(3,i), k*r0(4,i)
      end do

      return
      end
